Characterization and Demonstration of Mock Communities as Control Reagents for Accurate Human Microbiome Community Measurements

ABSTRACT Standardization and quality assurance of microbiome community analysis by high-throughput DNA sequencing require widely accessible and well-characterized reference materials. Here, we report on newly developed DNA and whole-cell mock communities to serve as control reagents for human gut microbiota measurements by shotgun metagenomics and 16S rRNA gene amplicon sequencing. The mock communities were formulated as near-even blends of up to 20 bacterial species prevalent in the human gut, span a wide range of genomic guanine-cytosine (GC) contents, and include multiple strains with Gram-positive type cell walls. Through a collaborative study, we carefully characterized the mock communities by shotgun metagenomics, using previously developed standardized protocols for DNA extraction and sequencing library construction. Further, we validated fitness of the mock communities for revealing technically meaningful differences among protocols for DNA extraction and metagenome/16S rRNA gene amplicon library construction. Finally, we used the mock communities to reveal varying performance of metagenome-based taxonomic profilers and the impact of trimming and filtering of sequencing reads on observed species profiles. The latter showed that aggressive preprocessing of reads may result in substantial GC-dependent bias and should thus be carefully evaluated to minimize unintended effects on species abundances. Taken together, the mock communities are expected to support a myriad of applications that rely on well-characterized control reagents, ranging from evaluation and optimization of methods to assessment of reproducibility in interlaboratory studies and routine quality control. IMPORTANCE Application of high-throughput DNA sequencing has greatly accelerated human microbiome research and its translation into new therapeutic and diagnostic capabilities. Microbiome community analyses results can, however, vary considerably across studies or laboratories, and establishment of measurement standards to improve accuracy and reproducibility has become a priority. The here-developed mock communities, which are available from the NITE Biological Resource Center (NBRC) at the National Institute of Technology and Evaluation (NITE, Japan), provide well-characterized control reagents that allow users to judge the accuracy of their measurement results. Widespread and consistent adoption of the mock communities will improve reproducibility and comparability of microbiome community analyses, thereby supporting and accelerating human microbiome research and development.


Protocols (non-SOPs) for DNA extraction
For Laboratory C, extraction of DNA was performed using the DNeasy PowerSoil HTP 96 Kit (Qiagen), following the manufacturer's recommended protocol. Bead-beating was performed using a MM 400 Mixer Mill (Retsch GmbH) at approximately 1,600 rpm for 20 min.
For Laboratory D, extraction of DNA was performed by a phenol/bead-beating method. First, 300 mg of 0.1mm glass beads, 300 μl of extraction buffer (167 mM Tris-HCl, pH 9.0, 67 mM EDTA and 1.7% SDS) and 500 μl of Tris-EDTA saturated phenol (Nippon Gene Co. Ltd., Tokyo, Japan) were added to the sample. Cells were then mechanically disrupted by bead-beating for 30 s at a power level of 5.0 in a FastPrep 24 instrument (MP Biomedicals, Irvine, CA, USA). After centrifugation at 15,000 × g for 5 min at 4ºC, the supernatant was transferred to a new 2.0-ml screw cap tube and mixed with 400 μl of phenol/chloroform/isoamyl alcohol (Nippon Gene Co. Ltd., Tokyo, Japan). The mixture was vortexed for 45 s at a power level of 4.0 using a FastPrep 24 and centrifuged at 15,000 × g for 5 min at 4ºC. The supernatant was then collected and DNA was precipitated by the addition 25 μl of 3 M sodium acetate (Nippon Gene Co. Ltd., Tokyo, Japan) and 300 μl of 2propanol (Sigma-Aldrich Co. LLC., St. Louis, MO, USA). After incubation at -30°C for approximately 60 min, the mixture was centrifuged at 15,000 × g for 5 min at 4ºC, and the pellet washed with 500 μl of 70% ethanol, and dissolved in 200 μl of Tris-EDTA buffer.
For Laboratory E, extraction of DNA was performed by combining 300 μl of lysis buffer (No. 10, Kurabo Industries Ltd.) and 0.5 g of 0.1-mm glass beads were mixed with the sample, followed by mechanical disruption using a Cell Destroyer PS1000 (Bio Medical Science, Tokyo, Japan) at approximately 4,000 rpm for 50 s at room temperature. Supernatant was then collected after centrifugation at 13,000 × g for 5 min and treated with 150 μl of proteinase K buffer (No. 2, Kurabo Industries Ltd; containing 0.4 mg ml -1 of proteinase K).
Purification of DNA was then performed with a Gene Prep Star PI-80X system (Kurabo Industries Ltd), following the manufacturer's provided instructions.

Protocols (non-SOPs) for metagenome library construction
For Laboratory B, shotgun metagenome libraries were constructed using the TruSeq DNA PCR-Free Kit (Illumina). To this end, DNA (1 μg) was fragmented by focused ultrasonication using a Covaris S220 instrument, with a target insert size of 550 bp. Fragmented DNA was then cleaned up and subjected to end-repair, adapter ligation and size selection, 3'-end adenylation and adapter ligation, following provided procedures.
For Laboratory C, shotgun metagenome libraries were prepared using an Illumina DNA Prep Kit (formerly, Nextera XT DNA Library Prep Kit), starting from 200 pg of input DNA, following manufacturer's provided instructions.
Thermal cycling conditions were as follows: 94ºC for 2 min; 18 cycles of 98ºC for 10 s, 55ºC for 30 s and 68ºC for 30 s. Triplicate PCR reactions were then pooled and purified using SPRIselect magnetic beads (Beckman Coulter). Amplicons were eluted with Buffer EB (Qiagen) and quantified by Quant-iT PicoGreen dsDNA Reagent (Thermo Fisher). To attach dual indexes and sequencing adapters, second-round PCR reactions (50 µl) using KOD -Plus-Ver.2 with 180 nM each of forward and reverse primers and 5 ng of purified first-round PCR products. Thermal cycling conditions were as follows: 94ºC for 2 min; 8 cycles of 98ºC for 10 s, 55ºC for 30 s and 68ºC for 45 s. Amplicons were then purified using SPRIselect magnetic beads and eluted with Buffer EB.
The DNA concentration was quantified by Quant-iT PicoGreen dsDNA Reagent.
The Nextera XT Index Kit was then used to attach dual indexes and sequencing adapters, in PCR reactions (20 µl) containing 1× Ex Taq buffer (Takara Bio), 200 μM of each dNTP (Takara Bio), 2 µl each of Index 1 and 2 primers, 0.5 units of TaKaRa ExTaq DNA polymerase and 0.2 μl of template DNA. Thermal cycling conditions were as follows: 98ºC for 30 s; 12 cycles of 98ºC for 40 s, 65ºC for 30 s and 72ºC for 30 s; 72ºC for 5 min.
Amplicons were then separated on 1% agarose gel and purified using the MinElute Gel Extraction Kit (Qiagen).
For Laboratory D, first round PCR reactions (25 μl) contained 0.625 units of Takara Ex Taq HS polymerase (Takara Bio, Shiga, Japan), 1 × Ex Taq Buffer, 0.2 mM each of dNTP, 200 nM each of forward (5'-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTGYCAGCMGCCGCGGTAA-3') and reverse primer (5'-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGACTACNVGGGTWTCTAAT-3') and 4.0 ng template DNA. Thermal cycling conditions were as follows: 95ºC for 3 min; 25 cycles of 95ºC for 30 s, 55ºC for 30 s and 72ºC for 30 s; 72ºC for 5 min. The amplified products were then used as templates for second round PCR to attach dual indexes and sequencing adapters, in PCR reactions (50 µl) containing 0.25 units of Takara Ex Taq HS polymerase (Takara Bio, Shiga, Japan), 1 × Ex Taq Buffer, 0.2 mM each of dNTP, 5.0 μl each of forward and reverse indexed primers of Nextera XT Index kit (Illumina, San Diego, CA, USA), and 1.0 µl of first-round PCR products. Thermal cycling conditions were as follows: 95ºC for 3 min; 10 cycles of 95ºC for 30 s, 55ºC for 30 s and 72ºC for 30 s; 72ºC for 5 min. Amplicons were purified using the Agencourt AMPure XP PCR Purification system and eluted in 10 mM Tris-Cl (pH8.5), and quantified using a Quant-iT PicoGreen dsDNA Assay kit (Life Technologies, Carlsbad, CA, USA).
Amplicon were then purified using AMPure XP (Beckman Coulter, Inc.) and eluted in 50 μl of 10 mM Tris-HCl, pH 8.5. The Nextera XT Index Kit Set A was then used to attach dual indexes and sequencing adapters, in PCR reactions (50 µl) containing 1× PCR buffer, 0.2 mM of each dNTP, 1 mM MgSO4, 5 µl each of forward and reverse indexed primers, 1 unit of KOD-Plus-v2 DNA polymerase (Toyobo Co., Ltd.) and 5 μl of template DNA.
Thermal cycling conditions were as follows: 95ºC for 3 min; 8 cycles of 95ºC for 30 s, 55ºC for 30 s and 68ºC for 1 min; and 68ºC for 5 min. Amplicons were then purified using the AMPure XP (Beckman Coulter, Inc.) and quantified with QuantiFluor ONE dsDNA System (Promega Co.).            Table S3. Figure S2. Taxonomic profiles of individual replicates for the DNA mock community (n = 16, upper facet) and cell mock community (n = 20, lower facet), generated by MetaPhlAn3. Shotgun metagenome libraries were prepared by two or three laboratories following SOPs for DNA extraction (protocol N) and sequencing library construction (protocols B and K), as indicated in the sample keys on the y axis. Fill colors and corresponding values overlayed onto the heatmap reflect relative abundances, as percentages.    Figure S4. Application of guidelines values for assessing methods for libray construction (DNA mock community, top panel) and DNA extraction + libray construction (cell mock community, lower panel). Blue and orange symbols represent the geometric mean and maximum of strain-wise absolute fold differences, respectively; corresponding "acceptance" values (see Table S4) are shown as horizontal lines with the same colors. Letters indicate participating laboratories. Numbers plotted at the top show the number of strains satisfying the "acceptance" range for strain-wise abundances (see Fig.1C).

A A A A A A A A A A A A B B B B B B B B B B B B C C C C C C C C C C C C
18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 17 17 17 17  Note that only the top-1000 most abundant species (based on the geometric mean) are shown. Data are shown as the geometric mean (circles) and standard deviation (error bars) of 16 measurements for the DNA mock community, performed by three laboratories following SOPs for DNA extraction and sequencing library construction. Grey curved segments connect species belonging to the same genus, originating from the most abundant species for a given genus.  Figure S6, cont'd. Same for kraken2_refseq. Figure S7. Agreement between the expected species profiles (relative abundances) and profiles generated by different taxonomic profilers. Agreement was calculated as Aitchison distances or Bray-Curtis dissimilarities between measured and expected species profiles. Note that only subcompositions of expected species (based on string matches to the assigned taxonomy strings) were considered, with rescaling to 100% for calculation of Bray-Curtis dissimilarities and replacement of zeros with 0.1% for calculation of Aitchison distances. Data are calculated for 16 individual measurements of the DNA mock community as performed using the SOPs (shapes) in three laboratories (colors). For the boxplots, the tick horizontal line represents the median, hinges show the 25th and 75th percentiles, whiskers extend to the largest and smallest value at most 1.5× the IQR (interquartile range) from the upper and lower hinges, respectively. Individual datapoints are overlayed with jitter.  Figure S8. Bray-Curtis dissimilarities between randomly subsampled datasets as generated using different profilers. To this end, we randomly subsampled four deeply sequenced HiSeq libraries, for the DNA mock community generated using protocols B and K (SRR15195689, SRR15195686, SRR15195688 and SRR15195685 in Table S3), to varying sequencing depths (0.5 to 20 million read pairs per sample). Subsequently, for each library and subsampling depth, all possible (45 per library, or 180 total) pairwise distances/dissimilarities between 10 random subsamples for each subsampling depth, as plotted on the x axis. For Bray-Curtis dissimilaties, all identified species, including false positives, were considered. For Aitchison distances, only expected species were considered.  Table S3, Q30 bases of 90.1%), medium (middle; library metagenome_mockCell_labB_lib043, Q30 bases of 78.2%) and low (right; library metagenome_mockCell_labB_lib039, Q30 bases of 71.7%) raw base quality are shown. Values shown on the y axis represent the ratio of retained reads for a given fastp setting (x axis) to the counts observed for the most lenient settings for trimming and filtering (that is, settings 0_100_50 in Table S5). Each line represents a different strain, with colors indicating the genomic GC content (%) as indicated in the legend.  Figure S10. Effect of read trimming and filtering using Cutadapt on strain-wise abundances. To this end, we subjected reads from three representative libraries with varying raw base call to initial trimming using fastp (settings 0_100_50 in Table S5) and then performed further quality trimming using Cutadapt for a range of quality thresholds as depicted on the x axis. Libraries are as in Figure S9, namely library metagenome_mockCell_labA_lib038 (Q30 bases of 90.1%, shown in blue), library metagenome_mockCell_labB_lib043 (Q30 bases of 78.2%, shown in orange) and library metagenome_mockCell_labB_lib039 (Q30 bases of 71.7%, shown in red). Facets (strains) are sorted by increasing genomic GC content.  Figure S11. Effect of quality trimming on number of classified species. For each profiler, the number of classified species was calculated for decreasing abundance thresholds as in Figure 3A for two different fastp settings. Plots represent scatter plots of the number of classified species for each abundance threshold for datasets with quality trimming (x axis) and without quality trimming (y axis). Lines are colored according to the percentage of bases with quality scores of >30 as indicated in the legend.  Figure S12. Distribution of strain-wise expected errors (EE), as calculated using USERACH's fastq_filter command, for the amplicon sequencing data; all available data for the DNA mock community (see Table S3) were included. For the boxplots, the tick horizontal line represents the median, hinges show the 25th and 75th percentiles, whiskers extend to the largest and smallest value at most 1.5× the IQR (interquartile range) from the upper and lower hinges, respectively; outliers are not shown. Note that participants A through D employed V2 chemistry for sequencing on a MiSeq instrument and participant E used V3 chemistry.  Figure S13. Raw base quality scores for sequences assigned to P. distasonis NBRC 113806. Data are shown as the mean (solid lines) and standard deviation (ribbons, if visible) of all available data for the DNA mock community (see Table S3). Note that participants A through D employed V2 chemistry for sequencing on a MiSeq instrument and participant E used V3 chemistry.